# 04_PSM_Rating.R
# Created: 2021-05-03 Taka-aki Asano
# Last Modified: 2024-02-16

# package
require("readr")
require("dplyr")
require("tidyr")
require("ggplot2")
require("ggpubr")
require("cjoint")
require("cregg")
theme_set(theme_minimal(base_size = 12))


# load dataset
data <- read_csv("HLWconjoint20210330.csv")


# processing of data
data2 <- select(data, 
                no, time, rl, group, age, psm01, psm02, psm03, psm04, 
                psm05, psm06, psm07, psm08, psm09, psm10, psm11, 
                psm12, psm13, psm14, psm15, psm16, hlw_su, hlw_fsu, 
                hlw_dis, hlw_pol, hlw_so, hlw_ad, hlw_pop, hlw_inc)
data2$hlw_fsu <- ifelse(data2$hlw_fsu == data2$rl, 1, 0)
data2 <- mutate(data2, 
                `Distance` = recode(hlw_dis, `1` = "less than 1 km", 
                                    `2` = "1-5 km", `3` = "5-20 km", 
                                    `4` = "20 km or more"), 
                `Politics` = recode(hlw_pol, `1` = "job creation", 
                                    `2` = "public services", 
                                    `3` = "local economy", 
                                    `4` = "society"), 
                `Social` = recode(hlw_so, `1` = "70% agreed, 30% disagreed", 
                                  `2` = "50% agreed, 50% disagreed", 
                                  `3` = "30% agreed, 70% disagreed"), 
                `Administration` = recode(hlw_ad, `1` = "none", 
                                          `2` = "childbirth and childcare", 
                                          `3` = "elderly care", `4` = "education", 
                                          `5` = "public facility", `6` = "water charge"), 
                `Population` = recode(hlw_pop, `1` = "19,000", 
                                      `2` = "20,000", `3` = "21,000"), 
                `Income` = recode(hlw_inc, `1` = "3.8 million yen", 
                                  `2` = "4 million yen", `3` = "4.2 million yen"))
data2$Distance <- factor(
  data2$Distance, 
  levels = c("less than 1 km", "1-5 km", "5-20 km", "20 km or more")
)
data2$Politics <- factor(
  data2$Politics, 
  levels = c("society", "job creation", "public services", "local economy")
)
data2$Social <- factor(
  data2$Social, 
  levels = c("30% agreed, 70% disagreed", "50% agreed, 50% disagreed", "70% agreed, 30% disagreed")
)
data2$Administration <- factor(
  data2$Administration, 
  levels = c("none", "childbirth and childcare", "elderly care", 
             "education", "public facility", "water charge")
)
data2$Population <- factor(
  data2$Population, 
  levels = c("19,000", "20,000", "21,000")
)
data2$Income <- factor(
  data2$Income, 
  levels = c("3.8 million yen", "4 million yen", "4.2 million yen")
)


# group
## PSM
data2$PSM <- rowSums(data2[,paste0("psm", c("01", "02", "03", "04", "05", "06", 
                                            "07", "08", "09", "10", "11", "12", 
                                            "13", "14", "15", "16"))])
summary(data2$PSM)
data2$PSM_Group <- ifelse(data2$PSM >= median(data2$PSM), "HPSM", "LPSM")
data2$PSM_Group <- factor(data2$PSM_Group, levels = c("LPSM", "HPSM"))
aggregate(PSM ~ PSM_Group, data2, mean)

## APS
data2$APS <- rowSums(data2[,paste0("psm", c("01", "03", "12", "15"))])
summary(data2$APS)
data2$APS_Group <- ifelse(data2$APS >= median(data2$APS), "HAPS", "LAPS")
data2$APS_Group <- factor(data2$APS_Group, levels = c("LAPS", "HAPS"))
aggregate(APS ~ APS_Group, data2, mean)

## CPV
data2$CPV <- rowSums(data2[,paste0("psm", c("04", "06", "07", "08"))])
summary(data2$CPV)
data2$CPV_Group <- ifelse(data2$CPV >= median(data2$CPV), "HCPV", "LCPV")
data2$CPV_Group <- factor(data2$CPV_Group, levels = c("LCPV", "HCPV"))
aggregate(CPV ~ CPV_Group, data2, mean)

## COM
data2$COM <- rowSums(data2[,paste0("psm", c("02", "09", "10", "16"))])
summary(data2$COM)
data2$COM_Group <- ifelse(data2$COM >= median(data2$COM), "HCOM", "LCOM")
data2$COM_Group <- factor(data2$COM_Group, levels = c("LCOM", "HCOM"))
aggregate(COM ~ COM_Group, data2, mean)

## SS
data2$SS <- rowSums(data2[,paste0("psm", c("05", "11", "13", "14"))])
summary(data2$SS)
data2$SS_Group <- ifelse(data2$SS >= median(data2$SS), "HSS", "LSS")
data2$SS_Group <- factor(data2$SS_Group, levels = c("LSS", "HSS"))
aggregate(SS ~ SS_Group, data2, mean)


# PSM
## AMCEs
res1 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ PSM_Group, level_order = "descending"
)
res1$level <- as.character(res1$level)
res1 <- rbind(
  res1, 
  c("HPSM", "hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HPSM")
)
res1$level <- factor(
  res1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res1$estimate <- as.numeric(res1$estimate)
res1$std.error <- as.numeric(res1$std.error)
res1$z <- as.numeric(res1$z)
res1$p <- as.numeric(res1$p)
res1$lower <- as.numeric(res1$lower)
res1$upper <- as.numeric(res1$upper)
res1.plot <- ggplot(
  res1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = PSM_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.55, 0.55) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res1.plot)

diff1 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ PSM_Group, level_order = "descending"
)
diff1$level <- as.character(diff1$level)
diff1 <- rbind(
  diff1, 
  c("HPSM - LPSM", "hlw_su", "amce", "Distance", "(Distance attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_su", "amce", "Politics", "(Political attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_su", "amce", "Social", "(Social attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_su", "amce", "Administration", "(Administrative attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_su", "amce", "Population", "(Population attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_su", "amce", "Income", "(Average income attribute)", "HPSM", NA, NA, NA, NA, NA, NA)
)
diff1$level <- factor(
  diff1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff1$estimate <- as.numeric(diff1$estimate)
diff1$std.error <- as.numeric(diff1$std.error)
diff1$z <- as.numeric(diff1$z)
diff1$p <- as.numeric(diff1$p)
diff1$lower <- as.numeric(diff1$lower)
diff1$upper <- as.numeric(diff1$upper)
diff1.plot <- ggplot(
  diff1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.5, 0.5) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff1.plot)


# APS
## AMCEs
res2 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ APS_Group, level_order = "descending"
)
res2$level <- as.character(res2$level)
res2 <- rbind(
  res2, 
  c("HAPS", "hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HAPS")
)
res2$level <- factor(
  res2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res2$estimate <- as.numeric(res2$estimate)
res2$std.error <- as.numeric(res2$std.error)
res2$z <- as.numeric(res2$z)
res2$p <- as.numeric(res2$p)
res2$lower <- as.numeric(res2$lower)
res2$upper <- as.numeric(res2$upper)
res2.plot <- ggplot(
  res2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = APS_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.55, 0.55) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res2.plot)

diff2 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ APS_Group, level_order = "descending"
)
diff2$level <- as.character(diff2$level)
diff2 <- rbind(
  diff2, 
  c("HAPS - LAPS", "hlw_su", "amce", "Distance", "(Distance attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_su", "amce", "Politics", "(Political attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_su", "amce", "Social", "(Social attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_su", "amce", "Administration", "(Administrative attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_su", "amce", "Population", "(Population attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_su", "amce", "Income", "(Average income attribute)", "HAPS", NA, NA, NA, NA, NA, NA)
)
diff2$level <- factor(
  diff2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff2$estimate <- as.numeric(diff2$estimate)
diff2$std.error <- as.numeric(diff2$std.error)
diff2$z <- as.numeric(diff2$z)
diff2$p <- as.numeric(diff2$p)
diff2$lower <- as.numeric(diff2$lower)
diff2$upper <- as.numeric(diff2$upper)
diff2.plot <- ggplot(
  diff2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.5, 0.5) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff2.plot)


# CPV
## AMCEs
res3 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ CPV_Group, level_order = "descending"
)
res3$level <- as.character(res3$level)
res3 <- rbind(
  res3, 
  c("HCPV", "hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCPV")
)
res3$level <- factor(
  res3$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res3$estimate <- as.numeric(res3$estimate)
res3$std.error <- as.numeric(res3$std.error)
res3$z <- as.numeric(res3$z)
res3$p <- as.numeric(res3$p)
res3$lower <- as.numeric(res3$lower)
res3$upper <- as.numeric(res3$upper)
res3.plot <- ggplot(
  res3, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = CPV_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.55, 0.55) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res3.plot)

diff3 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ CPV_Group, level_order = "descending"
)
diff3$level <- as.character(diff3$level)
diff3 <- rbind(
  diff3, 
  c("HCPV - LCPV", "hlw_su", "amce", "Distance", "(Distance attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_su", "amce", "Politics", "(Political attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_su", "amce", "Social", "(Social attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_su", "amce", "Administration", "(Administrative attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_su", "amce", "Population", "(Population attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_su", "amce", "Income", "(Average income attribute)", "HCPV", NA, NA, NA, NA, NA, NA)
)
diff3$level <- factor(
  diff3$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff3$estimate <- as.numeric(diff3$estimate)
diff3$std.error <- as.numeric(diff3$std.error)
diff3$z <- as.numeric(diff3$z)
diff3$p <- as.numeric(diff3$p)
diff3$lower <- as.numeric(diff3$lower)
diff3$upper <- as.numeric(diff3$upper)
diff3.plot <- ggplot(
  diff3, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.5, 0.5) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff3.plot)


# COM
## AMCEs
res4 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ COM_Group, level_order = "descending"
)
res4$level <- as.character(res4$level)
res4 <- rbind(
  res4, 
  c("HCOM", "hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCOM")
)
res4$level <- factor(
  res4$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res4$estimate <- as.numeric(res4$estimate)
res4$std.error <- as.numeric(res4$std.error)
res4$z <- as.numeric(res4$z)
res4$p <- as.numeric(res4$p)
res4$lower <- as.numeric(res4$lower)
res4$upper <- as.numeric(res4$upper)
res4.plot <- ggplot(
  res4, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = COM_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.55, 0.55) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res4.plot)

diff4 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ COM_Group, level_order = "descending"
)
diff4$level <- as.character(diff4$level)
diff4 <- rbind(
  diff4, 
  c("HCOM - LCOM", "hlw_su", "amce", "Distance", "(Distance attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_su", "amce", "Politics", "(Political attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_su", "amce", "Social", "(Social attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_su", "amce", "Administration", "(Administrative attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_su", "amce", "Population", "(Population attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_su", "amce", "Income", "(Average income attribute)", "HCOM", NA, NA, NA, NA, NA, NA)
)
diff4$level <- factor(
  diff4$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff4$estimate <- as.numeric(diff4$estimate)
diff4$std.error <- as.numeric(diff4$std.error)
diff4$z <- as.numeric(diff4$z)
diff4$p <- as.numeric(diff4$p)
diff4$lower <- as.numeric(diff4$lower)
diff4$upper <- as.numeric(diff4$upper)
diff4.plot <- ggplot(
  diff4, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.5, 0.5) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff4.plot)


# SS
## AMCEs
res5 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ SS_Group, level_order = "descending"
)
res5$level <- as.character(res5$level)
res5 <- rbind(
  res5, 
  c("HSS", "hlw_su", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_su", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_su", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_su", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_su", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_su", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HSS")
)
res5$level <- factor(
  res5$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res5$estimate <- as.numeric(res5$estimate)
res5$std.error <- as.numeric(res5$std.error)
res5$z <- as.numeric(res5$z)
res5$p <- as.numeric(res5$p)
res5$lower <- as.numeric(res5$lower)
res5$upper <- as.numeric(res5$upper)
res5.plot <- ggplot(
  res5, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = SS_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.55, 0.55) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res5.plot)

diff5 <- cj(
  data2, 
  hlw_su ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ SS_Group, level_order = "descending"
)
diff5$level <- as.character(diff5$level)
diff5 <- rbind(
  diff5, 
  c("HSS - LSS", "hlw_su", "amce", "Distance", "(Distance attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_su", "amce", "Politics", "(Political attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_su", "amce", "Social", "(Social attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_su", "amce", "Administration", "(Administrative attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_su", "amce", "Population", "(Population attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_su", "amce", "Income", "(Average income attribute)", "HSS", NA, NA, NA, NA, NA, NA)
)
diff5$level <- factor(
  diff5$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff5$estimate <- as.numeric(diff5$estimate)
diff5$std.error <- as.numeric(diff5$std.error)
diff5$z <- as.numeric(diff5$z)
diff5$p <- as.numeric(diff5$p)
diff5$lower <- as.numeric(diff5$lower)
diff5$upper <- as.numeric(diff5$upper)
diff5.plot <- ggplot(
  diff5, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.5, 0.5) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff5.plot)
